R Notebook

The goal here is to use conventional alpha diversity metrics to see how Chao1 richness, shannon diversity and evenness change across samples and to compare those to the values seen using breakaway in the AlphaDiversity.Rmd file

Setup

Run AlphaDiversity in scratchnotebooks That file calculates richness in breakawy which I will combine here

#source(here::here("RScripts", "InitialProcessing_3.R"))
source(here::here("RLibraries", "ChesapeakePersonalLibrary.R"))
ksource(here::here("ActiveNotebooks", "BreakawayAlphaDiversity.Rmd"))

processing file: /home/jacob/Projects/ChesapeakeMainstemAnalysis_ToShare/ActiveNotebooks/BreakawayAlphaDiversity.Rmd

  |                                                                                                           
  |                                                                                                     |   0%
  |                                                                                                           
  |....                                                                                                 |   4%
  |                                                                                                           
  |.........                                                                                            |   9%
  |                                                                                                           
  |.............                                                                                        |  13%
  |                                                                                                           
  |..................                                                                                   |  17%
  |                                                                                                           
  |......................                                                                               |  22%
  |                                                                                                           
  |..........................                                                                           |  26%
  |                                                                                                           
  |...............................                                                                      |  30%
  |                                                                                                           
  |...................................                                                                  |  35%
  |                                                                                                           
  |........................................                                                             |  39%
  |                                                                                                           
  |............................................                                                         |  43%
  |                                                                                                           
  |................................................                                                     |  48%
  |                                                                                                           
  |.....................................................                                                |  52%
  |                                                                                                           
  |.........................................................                                            |  57%
  |                                                                                                           
  |.............................................................                                        |  61%
  |                                                                                                           
  |..................................................................                                   |  65%
  |                                                                                                           
  |......................................................................                               |  70%
  |                                                                                                           
  |...........................................................................                          |  74%
  |                                                                                                           
  |...............................................................................                      |  78%
  |                                                                                                           
  |...................................................................................                  |  83%
  |                                                                                                           
  |........................................................................................             |  87%
  |                                                                                                           
  |............................................................................................         |  91%
  |                                                                                                           
  |.................................................................................................    |  96%
  |                                                                                                           
  |.....................................................................................................| 100%

output file: /tmp/RtmpswludR/file3e673faf9754

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

Attaching package: ‘flextable’

The following object is masked from ‘package:purrr’:

    compose


Attaching package: ‘ftExtra’

The following object is masked from ‘package:flextable’:

    separate_header

Warning: Assuming taxa are rows

library(vegan)

Loading required package: permute
Loading required package: lattice
This is vegan 2.6-3

library(cowplot)
library(flextable)
library(ftExtra)

This file is dedicated to conventional, non div-net/breakaway stats, since breakaway seems to choke on this data.

Reshape back into an ASV matrix, but this time correcting for total abundance

raDf <- nonSpikes_Remake %>% pivot_wider(id_cols = ID, names_from = ASV, values_from = RA, values_fill = 0)
raMat <- raDf %>% column_to_rownames("ID")

raMat1 <- raMat %>% as.matrix()

countMat <-  nonSpikes_Remake %>%
  pivot_wider(id_cols = ID, names_from = ASV, values_from = reads, values_fill = 0) %>%
  column_to_rownames("ID") %>% as.matrix()

min(seqDep)

[1] 852

This value is lower than the lowist chimera checked value because the spikes have been discarded (while chimera checked read depth still has spikes)

sampleRichness <- rarefy(countMat, min(seqDep))

rarefy everything to the minimum depth (852)

countRare <- rrarefy(countMat, min(seqDep))

Gamma diversity

specpool(countRare)

Doesn’t finish

#specpool(countMat)

Calculate diversity indeces

All richness estimates

richnessRare <- estimateR(countRare)

Shannon diversity

shan <- diversity(countRare)
shan

 3-1-B-0-2  3-1-B-1-2  3-1-B-180   3-1-B-20    3-1-B-5  3-1-B-500   3-1-B-53  3-1-S-0-2  3-1-S-1-2  3-1-S-180 
  4.447712   5.116637   4.678407   5.866088   5.097057   3.871649   5.479328   4.615634   4.929580   4.704532 
  3-1-S-20    3-1-S-5  3-2-B-0-2  3-2-B-1-2  3-2-B-180   3-2-B-20    3-2-B-5  3-2-B-500   3-2-B-53  3-2-S-0-2 
  5.260828   4.837697   4.502874   4.691918   4.620095   5.164912   5.416675   4.988461   4.938664   3.732937 
 3-2-S-1-2  3-2-S-180   3-2-S-20    3-2-S-5  3-2-S-500   3-2-S-53  3-3-B-0-2  3-3-B-1-2  3-3-B-180   3-3-B-20 
  4.849386   4.751104   4.843166   5.166633   4.808113   4.307320   4.350618   4.971685   3.437542   5.685777 
   3-3-B-5  3-3-B-500   3-3-B-53  3-3-S-180   3-3-S-20  3-3-S-500   3-3-S-53  4-3-B-0-2  4-3-B-1-2  4-3-B-180 
  5.260559   5.262516   5.343435   4.981628   4.942967   4.859932   4.357677   4.296520   5.024896   4.583716 
  4-3-B-20    4-3-B-5  4-3-B-500   4-3-B-53  4-3-O-1-2  4-3-O-180    4-3-O-5  4-3-O-500   4-3-O-53  4-3-S-0-2 
  4.394405   4.630198   4.426126   4.227551   4.968388   4.715856   5.195087   4.034619   4.583134   2.805419 
 4-3-S-180   4-3-S-20  4-3-S-500   4-3-S-53  5-1-S-1-2  5-1-S-180   5-1-S-20    5-1-S-5  5-1-S-500   5-1-S-53 
  4.563265   4.627756   4.805569   4.460398   4.284061   4.554498   4.210847   3.956005   4.066923   4.196692 
 5-5-B-0-2  5-5-B-180    5-5-B-5  5-5-B-500   5-5-B-53  5-5-S-180    5-5-S-5  5-5-S-500   5-5-S-53 C_5P1B_0P2 
  4.573474   5.192896   5.419692   5.098760   4.829526   4.316809   5.003296   4.843684   4.216447   4.377762 
C_5P1B_180 C_5P1B_1P2  C_5P1B_20 C_5P1B_500  C_5P1B_53 
  4.814288   4.856735   5.376345   4.676670   5.054771 

Evenness

pielouJ <- shan/richnessRare["S.chao1",]
pielouJ

  3-1-B-0-2   3-1-B-1-2   3-1-B-180    3-1-B-20     3-1-B-5   3-1-B-500    3-1-B-53   3-1-S-0-2   3-1-S-1-2 
0.012242431 0.007628407 0.010461256 0.002668322 0.006980389 0.042083146 0.007982186 0.008122541 0.007524956 
  3-1-S-180    3-1-S-20     3-1-S-5   3-2-B-0-2   3-2-B-1-2   3-2-B-180    3-2-B-20     3-2-B-5   3-2-B-500 
0.007157256 0.008075588 0.006058481 0.012802536 0.008688737 0.013388595 0.005490156 0.007066346 0.009079896 
   3-2-B-53   3-2-S-0-2   3-2-S-1-2   3-2-S-180    3-2-S-20     3-2-S-5   3-2-S-500    3-2-S-53   3-3-B-0-2 
0.010109855 0.026272434 0.007768416 0.009249390 0.009327234 0.007932156 0.006879380 0.014943315 0.011706652 
  3-3-B-1-2   3-3-B-180    3-3-B-20     3-3-B-5   3-3-B-500    3-3-B-53   3-3-S-180    3-3-S-20   3-3-S-500 
0.007745782 0.061937694 0.005071165 0.006867398 0.006900477 0.005323821 0.009684942 0.008590692 0.008665187 
   3-3-S-53   4-3-B-0-2   4-3-B-1-2   4-3-B-180    4-3-B-20     4-3-B-5   4-3-B-500    4-3-B-53   4-3-O-1-2 
0.009491849 0.009200256 0.007283966 0.010619160 0.008898709 0.009710793 0.010432976 0.008489059 0.008707836 
  4-3-O-180     4-3-O-5   4-3-O-500    4-3-O-53   4-3-S-0-2   4-3-S-180    4-3-S-20   4-3-S-500    4-3-S-53 
0.012332810 0.005953768 0.016115035 0.012607328 0.124685269 0.012810964 0.008107808 0.010446890 0.013387723 
  5-1-S-1-2   5-1-S-180    5-1-S-20     5-1-S-5   5-1-S-500    5-1-S-53   5-5-B-0-2   5-5-B-180     5-5-B-5 
0.008619198 0.012719017 0.014697546 0.015157109 0.015985286 0.015201275 0.007584534 0.007587199 0.006467270 
  5-5-B-500    5-5-B-53   5-5-S-180     5-5-S-5   5-5-S-500    5-5-S-53  C_5P1B_0P2  C_5P1B_180  C_5P1B_1P2 
0.009734807 0.009680381 0.014795877 0.009898759 0.008377902 0.016317124 0.008321046 0.007915634 0.005932011 
  C_5P1B_20  C_5P1B_500   C_5P1B_53 
0.003437984 0.008991717 0.006222376 

Combine diversity data

diversityData <- sampleData %>% 
  left_join(richnessRare %>% t() %>% as.data.frame() %>% rownames_to_column("ID"), by = "ID") %>%
  left_join(shan %>% enframe(name = "ID", value = "shannonH"), by = "ID") %>%
  left_join(pielouJ %>% enframe(name = "ID", value = "pielouJ"), by = "ID") %>%
  arrange(Size_Class)

Generate plots of diversity estimates

Parameters for all plots

plotSpecs <- list(
  facet_wrap(~Depth, ncol = 1) ,
  theme_bw(base_size = 16) ,
  geom_point(size = 4) ,
  geom_path(aes(color = as.factor(Station))) ,
  scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
  #scale_y_log10nice() ,
  scale_shape_manual(values = rep(21:25, 2)) ,
  scale_fill_viridis_d(option = "plasma") ,
  scale_color_viridis_d(option = "plasma") ,
  labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
  theme(legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = .5),
        axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)

Observed species counts, on rarefied data

plotObs <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = S.obs, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +ylab("Observed ASVs (Rarefied)")#+ scale_y_log10()
plotObs

Estemated Chao1 Richness

plotChao1 <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = S.chao1, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  geom_errorbar(aes(ymin = S.chao1 -2 * se.chao1, ymax = S.chao1 + 2* se.chao1), width = -.1) + 
  scale_y_log10() +
  ylab("Richness (Chao1)")
plotChao1

Shannon diversity

plotShan <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = shannonH, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  ylab("Diversity (Shannon H)") +
  lims(y = c(2.5, 6))
plotShan

Evenness

plotPielou <- diversityData %>%

#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +scale_y_log10() +ylab("Evenness (PielouJ)")
plotPielou

All plots together

plotAlpha <- plot_grid(plotObs, plotChao1, plotShan, plotPielou, nrow = 1, labels = LETTERS)
plotAlpha

ggsave(here::here("Figures", "ConventionalAlpha.png"), plotAlpha, width = 11, height = 4)

Do we see trends with lat and size?

Observed Species

Rarefied observed species numbers

obsMod <- lm(S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(obsMod)

Call:
lm(formula = S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + 
    I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

Residuals:
     Min       1Q   Median       3Q      Max 
-221.223  -35.304    4.193   45.989  200.602 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)          253678.231 131330.970   1.932 0.057517 .  
log(Size_Class)          31.672      7.957   3.980 0.000168 ***
I(log(Size_Class)^2)     -6.197      1.480  -4.186 8.24e-05 ***
lat                  -13212.822   6843.487  -1.931 0.057628 .  
I(lat^2)                172.124     89.089   1.932 0.057461 .  
depth                     4.111      3.151   1.305 0.196386    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 74.8 on 69 degrees of freedom
Multiple R-squared:  0.2634,    Adjusted R-squared:   0.21 
F-statistic: 4.935 on 5 and 69 DF,  p-value: 0.0006408

Richness

Rarified chao1 estimates

chao1Mod <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(chao1Mod)

Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + 
    I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-544.01 -152.34  -28.53  144.59 1401.31 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)   
(Intercept)           1.037e+06  5.008e+05   2.071  0.04211 * 
log(Size_Class)       8.827e+01  3.034e+01   2.909  0.00487 **
I(log(Size_Class)^2) -1.931e+01  5.645e+00  -3.421  0.00105 **
lat                  -5.406e+04  2.610e+04  -2.072  0.04204 * 
I(lat^2)              7.044e+02  3.397e+02   2.074  0.04185 * 
depth                 2.076e+01  1.202e+01   1.727  0.08857 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 285.2 on 69 degrees of freedom
Multiple R-squared:  0.2082,    Adjusted R-squared:  0.1509 
F-statistic:  3.63 on 5 and 69 DF,  p-value: 0.005664

As above but without latitude and depth

chao1ModSimple <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), data = diversityData)
summary(chao1ModSimple)

Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), 
    data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-469.09 -166.48  -32.74  134.06 1512.72 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           597.962     51.919  11.517  < 2e-16 ***
log(Size_Class)        88.481     30.833   2.870 0.005391 ** 
I(log(Size_Class)^2)  -19.759      5.736  -3.445 0.000956 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 290.2 on 72 degrees of freedom
Multiple R-squared:  0.1447,    Adjusted R-squared:  0.121 
F-statistic: 6.092 on 2 and 72 DF,  p-value: 0.003596

Shannon Diversity

shanMod <- lm(shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

summary(shanMod)

Call:
lm(formula = shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityData)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.33049 -0.22559  0.06304  0.28864  0.70832 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.424e+03  7.757e+02   1.836 0.070648 .  
log(Size_Class)       1.914e-01  4.700e-02   4.072 0.000122 ***
I(log(Size_Class)^2) -3.501e-02  8.744e-03  -4.004 0.000155 ***
lat                  -7.395e+01  4.042e+01  -1.830 0.071627 .  
I(lat^2)              9.627e-01  5.262e-01   1.830 0.071630 .  
depth                 1.458e-02  1.861e-02   0.783 0.436226    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4418 on 69 degrees of freedom
Multiple R-squared:  0.2814,    Adjusted R-squared:  0.2294 
F-statistic: 5.405 on 5 and 69 DF,  p-value: 0.0002982

Evenness

pielouMod <- lm(pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(pielouMod)

Call:
lm(formula = pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityData)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.016051 -0.004926 -0.002385  0.001354  0.099434 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -2.747e+01  2.640e+01  -1.041   0.3016  
log(Size_Class)      -4.045e-03  1.599e-03  -2.529   0.0137 *
I(log(Size_Class)^2)  7.144e-04  2.976e-04   2.401   0.0191 *
lat                   1.431e+00  1.375e+00   1.040   0.3019  
I(lat^2)             -1.861e-02  1.791e-02  -1.039   0.3023  
depth                -4.148e-04  6.334e-04  -0.655   0.5148  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01504 on 69 degrees of freedom
Multiple R-squared:  0.1054,    Adjusted R-squared:  0.0406 
F-statistic: 1.626 on 5 and 69 DF,  p-value: 0.1646

uomisto H (2010a). “A diversity of beta diver- sities: straightening up a concept gone awry. 1. Defining beta diversity as a function of alpha and gamma diversity.” Ecography, 33, 2–2

Prediction plots

Observed Species

predict(obsMod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12 
230.1017 230.1017 201.2125 201.2125 209.8268 162.3905 162.3905 187.5034 200.5021 302.6958 302.6958 273.8066 
      13       14       15       16       17       18       19       20       21       22       23       24 
273.8066 282.4209 234.9846 234.9846 260.0975 260.0975 332.0508 332.0508 303.1616 303.1616 311.7759 264.3396 
      25       26       27       28       29       30       31       32       33       34       35       36 
264.3396 289.4525 302.4512 302.4512 336.3980 336.3980 307.5088 307.5088 316.1232 316.1232 268.6868 268.6868 
      37       38       39       40       41       42       43       44       45       46       47       48 
293.7998 293.7998 325.1970 296.3078 296.3078 304.9221 304.9221 257.4858 257.4858 257.4858 282.5987 282.5987 
      49       50       51       52       53       54       55       56       57       58       59       60 
295.5974 295.5974 294.4974 294.4974 265.6083 265.6083 274.2226 274.2226 226.7862 226.7862 226.7862 251.8992 
      61       62       63       64       65       66       67       68       69       70       71       72 
251.8992 264.8979 264.8979 254.6368 225.7477 225.7477 234.3620 234.3620 186.9256 186.9256 186.9256 212.0386 
      73       74       75 
212.0386 225.0373 225.0373 

$se.fit
 [1] 26.16949 26.16949 25.45599 25.45599 25.32015 26.93768 26.93768 28.45024 33.00461 18.35754 18.35754
[12] 17.75988 17.75988 16.67860 19.43184 19.43184 20.85030 20.85030 18.66951 18.66951 18.22171 18.22171
[23] 16.71287 19.47713 19.47713 20.56141 27.32163 27.32163 18.66740 18.66740 18.18962 18.18962 16.44692
[34] 16.44692 19.01561 19.01561 19.95080 19.95080 17.99853 17.37889 17.37889 15.51257 15.51257 17.88327
[45] 17.88327 17.88327 18.84420 18.84420 25.74850 25.74850 18.72596 18.72596 17.86591 17.86591 16.17413
[56] 16.17413 17.86737 17.86737 17.86737 18.90703 18.90703 25.32224 25.32224 23.71732 22.79314 22.79314
[67] 21.66683 21.66683 22.43854 22.43854 22.43854 23.40764 23.40764 28.33970 28.33970

$df
[1] 69

$residual.scale
[1] 74.80449

diversityData$pred_obs = predict(obsMod, se.fit = TRUE)$fit
diversityData$se_obs = predict(obsMod, se.fit = TRUE)$se.fit

plotSpecs2 <- list(
  facet_wrap(~Depth, ncol = 1) ,
  theme_bw(base_size = 16) ,
  #geom_point(size = 4) ,
  geom_path(aes(color = as.factor(Station))) ,
  scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
  #scale_y_log10nice() ,
  scale_shape_manual(values = rep(21:25, 2)) ,
  scale_fill_viridis_d(option = "plasma") ,
  scale_color_viridis_d(option = "plasma") ,
  labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
  theme(legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = .5),
        axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)

plotObs_pred <-  diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_obs, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_obs - 2 * se_obs, yend = pred_obs + 2 * se_obs, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5)  +
  plotSpecs2 + ylab("Predicted  ASVs") 
plotObs_pred

Richness

predict(chao1Mod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12 
513.9062 513.9062 390.1285 390.1285 469.8081 281.2073 281.2073 408.9487 387.9514 721.4337 721.4337 597.6560 
      13       14       15       16       17       18       19       20       21       22       23       24 
597.6560 677.3356 488.7348 488.7348 616.4761 616.4761 798.0222 798.0222 674.2445 674.2445 753.9241 565.3233 
      25       26       27       28       29       30       31       32       33       34       35       36 
565.3233 693.0646 672.0674 672.0674 797.1065 797.1065 673.3289 673.3289 753.0084 753.0084 564.4077 564.4077 
      37       38       39       40       41       42       43       44       45       46       47       48 
692.1490 692.1490 752.0344 628.2567 628.2567 707.9363 707.9363 519.3355 519.3355 519.3355 647.0769 647.0769 
      49       50       51       52       53       54       55       56       57       58       59       60 
626.0796 626.0796 643.6118 643.6118 519.8342 519.8342 599.5137 599.5137 410.9130 410.9130 410.9130 538.6543 
      61       62       63       64       65       66       67       68       69       70       71       72 
538.6543 517.6571 517.6571 508.7384 384.9608 384.9608 464.6403 464.6403 276.0396 276.0396 276.0396 403.7809 
      73       74       75 
403.7809 382.7836 382.7836 

$se.fit
 [1]  99.78989  99.78989  97.06920  97.06920  96.55117 102.71918 102.71918 108.48688 125.85368  70.00126
[11]  70.00126  67.72225  67.72225  63.59909  74.09781  74.09781  79.50669  79.50669  71.19088  71.19088
[21]  69.48331  69.48331  63.72976  74.27049  74.27049  78.40508 104.18328 104.18328  71.18284  71.18284
[31]  69.36094  69.36094  62.71564  62.71564  72.51061  72.51061  76.07670  76.07670  68.63229  66.26945
[41]  66.26945  59.15278  59.15278  68.19276  68.19276  68.19276  71.85699  71.85699  98.18459  98.18459
[51]  71.40613  71.40613  68.12657  68.12657  61.67546  61.67546  68.13214  68.13214  68.13214  72.09658
[61]  72.09658  96.55916  96.55916  90.43926  86.91517  86.91517  82.62028  82.62028  85.56300  85.56300
[71]  85.56300  89.25838  89.25838 108.06537 108.06537

$df
[1] 69

$residual.scale
[1] 285.2456

diversityData$pred_chao1 = predict(chao1Mod, se.fit = TRUE)$fit
diversityData$se_chao1 = predict(chao1Mod, se.fit = TRUE)$se.fit

plotChao1_pred <-  diversityData %>%

#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_chao1, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_chao1 - 2 * se_chao1, yend = pred_chao1 + 2 * se_chao1, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5)  +
  plotSpecs2 + ylab("Predictd Richness (Chao1)") + scale_y_log10()
plotChao1_pred

Shannon Diversity

predict(shanMod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12 
4.510911 4.510911 4.363021 4.363021 4.319694 4.051892 4.051892 4.143951 4.388759 4.943313 4.943313 4.795423 
      13       14       15       16       17       18       19       20       21       22       23       24 
4.795423 4.752096 4.484294 4.484294 4.576353 4.576353 5.126894 5.126894 4.979003 4.979003 4.935677 4.667874 
      25       26       27       28       29       30       31       32       33       34       35       36 
4.667874 4.759933 5.004741 5.004741 5.168679 5.168679 5.020788 5.020788 4.977462 4.977462 4.709659 4.709659 
      37       38       39       40       41       42       43       44       45       46       47       48 
4.801718 4.801718 5.117504 4.969614 4.969614 4.926287 4.926287 4.658485 4.658485 4.658485 4.750544 4.750544 
      49       50       51       52       53       54       55       56       57       58       59       60 
4.995352 4.995352 4.959250 4.959250 4.811359 4.811359 4.768033 4.768033 4.500230 4.500230 4.500230 4.592289 
      61       62       63       64       65       66       67       68       69       70       71       72 
4.592289 4.837097 4.837097 4.746740 4.598849 4.598849 4.555523 4.555523 4.287720 4.287720 4.287720 4.379779 
      73       74       75 
4.379779 4.624587 4.624587 

$se.fit
 [1] 0.15456512 0.15456512 0.15035102 0.15035102 0.14954865 0.15910231 0.15910231 0.16803594 0.19493546
[10] 0.10842534 0.10842534 0.10489537 0.10489537 0.09850898 0.11477052 0.11477052 0.12314836 0.12314836
[19] 0.11026796 0.11026796 0.10762308 0.10762308 0.09871139 0.11503798 0.11503798 0.12144207 0.16137006
[28] 0.16137006 0.11025550 0.11025550 0.10743355 0.10743355 0.09714060 0.09714060 0.11231209 0.11231209
[37] 0.11783563 0.11783563 0.10630494 0.10264512 0.10264512 0.09162207 0.09162207 0.10562415 0.10562415
[46] 0.10562415 0.11129969 0.11129969 0.15207866 0.15207866 0.11060136 0.11060136 0.10552162 0.10552162
[55] 0.09552947 0.09552947 0.10553025 0.10553025 0.10553025 0.11167080 0.11167080 0.14956102 0.14956102
[64] 0.14008188 0.13462339 0.13462339 0.12797101 0.12797101 0.13252901 0.13252901 0.13252901 0.13825280
[73] 0.13825280 0.16738306 0.16738306

$df
[1] 69

$residual.scale
[1] 0.4418186

diversityData$pred_shanH = predict(shanMod, se.fit = TRUE)$fit
diversityData$se_shanH = predict(shanMod, se.fit = TRUE)$se.fit

plotShannonH_pred <- diversityData %>%

 #filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_shanH, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_shanH - 2 * se_shanH, yend = pred_shanH + 2 * se_shanH, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"),  alpha = 0.5)  +
  plotSpecs2 + ylab("Predicted Diversity (Shannon H)") #+ scale_y_log10()
plotShannonH_pred

Evenness

predict(pielouMod, se.fit = TRUE)

$fit
          1           2           3           4           5           6           7           8           9 
0.019088298 0.019088298 0.021985268 0.021985268 0.021029079 0.025251377 0.025251377 0.022216580 0.018901598 
         10          11          12          13          14          15          16          17          18 
0.010013570 0.010013570 0.012910541 0.012910541 0.011954352 0.016176649 0.016176649 0.013141852 0.013141852 
         19          20          21          22          23          24          25          26          27 
0.006067283 0.006067283 0.008964253 0.008964253 0.008008064 0.012230362 0.012230362 0.009195565 0.005880583 
         28          29          30          31          32          33          34          35          36 
0.005880583 0.005020093 0.005020093 0.007917064 0.007917064 0.006960875 0.006960875 0.011183172 0.011183172 
         37          38          39          40          41          42          43          44          45 
0.008148375 0.008148375 0.005927552 0.008824522 0.008824522 0.007868333 0.007868333 0.012090631 0.012090631 
         46          47          48          49          50          51          52          53          54 
0.012090631 0.009055834 0.009055834 0.005740852 0.005740852 0.008985179 0.008985179 0.011882150 0.011882150 
         55          56          57          58          59          60          61          62          63 
0.010925961 0.010925961 0.015148258 0.015148258 0.015148258 0.012113461 0.012113461 0.008798480 0.008798480 
         64          65          66          67          68          69          70          71          72 
0.013178110 0.016075081 0.016075081 0.015118892 0.015118892 0.019341189 0.019341189 0.019341189 0.016306392 
         73          74          75 
0.016306392 0.012991411 0.012991411 

$se.fit
 [1] 0.005259890 0.005259890 0.005116483 0.005116483 0.005089178 0.005414292 0.005414292 0.005718306
 [9] 0.006633703 0.003689742 0.003689742 0.003569616 0.003569616 0.003352286 0.003905670 0.003905670
[17] 0.004190770 0.004190770 0.003752447 0.003752447 0.003662441 0.003662441 0.003359174 0.003914772
[25] 0.003914772 0.004132704 0.005491464 0.005491464 0.003752023 0.003752023 0.003655991 0.003655991
[33] 0.003305719 0.003305719 0.003822009 0.003822009 0.004009976 0.004009976 0.003617584 0.003493039
[41] 0.003493039 0.003117922 0.003117922 0.003594417 0.003594417 0.003594417 0.003787557 0.003787557
[49] 0.005175275 0.005175275 0.003763792 0.003763792 0.003590927 0.003590927 0.003250892 0.003250892
[57] 0.003591221 0.003591221 0.003591221 0.003800186 0.003800186 0.005089600 0.005089600 0.004767022
[65] 0.004581268 0.004581268 0.004354886 0.004354886 0.004509996 0.004509996 0.004509996 0.004704778
[73] 0.004704778 0.005696088 0.005696088

$df
[1] 69

$residual.scale
[1] 0.0150352

diversityData$pred_pielouJ = predict(pielouMod, se.fit = TRUE)$fit
diversityData$se_pielouJ = predict(pielouMod, se.fit = TRUE)$se.fit

plot_pielouJ_pred <- diversityData %>%

#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_pielouJ - 2 * se_pielouJ, yend = pred_pielouJ + 2 * se_pielouJ, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Evenness (Pielou J)") + scale_y_log10()
plot_pielouJ_pred

Combined prediction plot

plotPredictions <- plot_grid(plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred, nrow = 1, labels = LETTERS)

Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (`geom_segment()`).

plotPredictions

ggsave(here::here("Figures", "ConventionalAlphaPredictions.png"), plotPredictions, width = 11, height = 4)

Even combindeder

plot_grid(plotObs, plotChao1, plotShan, plotPielou,
          plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred,
          nrow = 2, labels = LETTERS)

Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (`geom_segment()`).

Combined summary table

alphaSummary <- tibble(
  metric = c("Observed ASVs", "Richness (Chao1)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
  model = list(obsMod, chao1Mod, shanMod, pielouMod)
)

alphaSummary <- alphaSummary %>%
  mutate(df = map(model, ~broom::tidy(summary(.))))

alphaSummary <- alphaSummary %>%
  select(-model) %>%
  unnest(df)

alphaSummary <- alphaSummary %>%
  rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
  mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
         Term = str_replace(Term, "\\^2", "\\^2\\^"),
         Term = str_replace(Term, "depth", "Depth"),
         Term = str_replace(Term, "lat", "Latitude"),
         Term = str_replace(Term, "_", " ")# BOOKMARK!!
         ) %>%
  mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
           reformat_sci()
         ) %>%
  mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
           reformat_sci()
  ) %>%
  mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
  mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
  rename(`Standard\nError` = `Standard Error`) %>%
  identity()

alphaSummary %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>%
  bold(i = ~ p< 0.05, j = "p") %>%
  colformat_md() %>%
  set_table_properties(layout = "autofit") %>%
  align(j = -c(1:2), align = "right")

Metric	Term	Estimate	Standard Error	T Value	p
Observed ASVs	Intercept	2.5 x 105	1.3 x 105	1.93	0.058
	log(Size Class)	3.2 x 101	8.0 x 100	3.98	< 0.001
	log(Size Class)2	-6.2 x 100	1.5 x 100	-4.19	< 0.001
	Latitude	-1.3 x 104	6.8 x 103	-1.93	0.058
	Latitude2	1.7 x 102	8.9 x 101	1.93	0.057
	Depth	4.1 x 100	3.2 x 100	1.30	0.196
Richness (Chao1)	Intercept	1.0 x 106	5.0 x 105	2.07	0.042
	log(Size Class)	8.8 x 101	3.0 x 101	2.91	0.005
	log(Size Class)2	-1.9 x 101	5.6 x 100	-3.42	0.001
	Latitude	-5.4 x 104	2.6 x 104	-2.07	0.042
	Latitude2	7.0 x 102	3.4 x 102	2.07	0.042
	Depth	2.1 x 101	1.2 x 101	1.73	0.089
Diversity (Shannon H)	Intercept	1.4 x 103	7.8 x 102	1.84	0.071
	log(Size Class)	1.9 x 10-1	4.7 x 10-2	4.07	< 0.001
	log(Size Class)2	-3.5 x 10-2	8.7 x 10-3	-4.00	< 0.001
	Latitude	-7.4 x 101	4.0 x 101	-1.83	0.072
	Latitude2	9.6 x 10-1	5.3 x 10-1	1.83	0.072
	Depth	1.5 x 10-2	1.9 x 10-2	0.78	0.436
Evenness (Pielou J)	Intercept	-2.7 x 101	2.6 x 101	-1.04	0.302
	log(Size Class)	-4.0 x 10-3	1.6 x 10-3	-2.53	0.014
	log(Size Class)2	7.1 x 10-4	3.0 x 10-4	2.40	0.019
	Latitude	1.4 x 100	1.4 x 100	1.04	0.302
	Latitude2	-1.9 x 10-2	1.8 x 10-2	-1.04	0.302
	Depth	-4.1 x 10-4	6.3 x 10-4	-0.65	0.515

Now considering breakaway values

richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.))))

Warning: `rename_()` was deprecated in dplyr 0.7.0.
Please use `rename()` instead.

diversityDataWB <- full_join(diversityData,
                             richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.)))),
                             by = c("ID" = "break_sample_names"), suffix = c("", "_break")) %>%
  mutate(pielouJ2 = shannonH/break_estimate) %>%
  identity()

diversityDataWB

pielouMod2 <- lm(pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityDataWB)
summary(pielouMod2)

Call:
lm(formula = pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityDataWB)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.013975 -0.005102 -0.002495  0.000889  0.105935 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -1.971e+01  2.752e+01  -0.716   0.4762  
log(Size_Class)      -3.292e-03  1.668e-03  -1.974   0.0524 .
I(log(Size_Class)^2)  5.752e-04  3.103e-04   1.854   0.0680 .
lat                   1.025e+00  1.434e+00   0.715   0.4771  
I(lat^2)             -1.332e-02  1.867e-02  -0.713   0.4780  
depth                -2.399e-04  6.605e-04  -0.363   0.7176  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01568 on 69 degrees of freedom
Multiple R-squared:  0.06743,   Adjusted R-squared:  -0.000147 
F-statistic: 0.9978 on 5 and 69 DF,  p-value: 0.4257

Ok. So the narrative makes sense. Alpha diveristy is driven by variability in richness rather than evenness. Why would we see an effect in chao1 but not breakaway? Because chao1 is more driven by abundant stuff that makes the rarification threshold. My first hunch is that chao1 responds to evenness, but actually that shouldn’t have any effect since there is no evenness variability? Or maybe just that breakaway is more variable (because it detects fine level differences in rare species) and that doesn’t map as nicely with overall patterns.

plotBreak <- diversityDataWB %>%

#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  #scale_y_log10()+
  ylab("Richness (Breakaway)")
plotBreak

plotPielou2 <- diversityDataWB %>%

#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  #scale_y_log10()+
  ylab("Evenness (PielouJ)")
plotPielou2

Redo predictions for good measure

predict(pielouMod2, se.fit = TRUE)

$fit
            1             2             3             4             5             6             7 
 0.0116808693  0.0116808693  0.0135879652  0.0135879652  0.0133789427  0.0160393203  0.0160393203 
            8             9            10            11            12            13            14 
 0.0139920245  0.0099168119  0.0043118512  0.0043118512  0.0062189471  0.0062189471  0.0060099245 
           15            16            17            18            19            20            21 
 0.0086703022  0.0086703022  0.0066230064  0.0066230064  0.0010848943  0.0010848943  0.0029919902 
           22            23            24            25            26            27            28 
 0.0029919902  0.0027829677  0.0054433453  0.0054433453  0.0033960495 -0.0006791631 -0.0006791631 
           29            30            31            32            33            34            35 
 0.0001937160  0.0001937160  0.0021008119  0.0021008119  0.0018917894  0.0018917894  0.0045521670 
           36            37            38            39            40            41            42 
 0.0045521670  0.0025048712  0.0025048712  0.0008906780  0.0027977739  0.0027977739  0.0025887513 
           43            44            45            46            47            48            49 
 0.0025887513  0.0052491290  0.0052491290  0.0052491290  0.0032018332  0.0032018332 -0.0008733794 
           50            51            52            53            54            55            56 
-0.0008733794  0.0033103787  0.0033103787  0.0052174746  0.0052174746  0.0050084520  0.0050084520 
           57            58            59            60            61            62            63 
 0.0076688297  0.0076688297  0.0076688297  0.0056215339  0.0056215339  0.0015463213  0.0015463213 
           64            65            66            67            68            69            70 
 0.0066511872  0.0085582831  0.0085582831  0.0083492605  0.0083492605  0.0110096382  0.0110096382 
           71            72            73            74            75 
 0.0110096382  0.0089623424  0.0089623424  0.0048871298  0.0048871298 

$se.fit
 [1] 0.005484417 0.005484417 0.005334888 0.005334888 0.005306418 0.005645409 0.005645409 0.005962400
 [9] 0.006916873 0.003847244 0.003847244 0.003721991 0.003721991 0.003495383 0.004072389 0.004072389
[17] 0.004369659 0.004369659 0.003912625 0.003912625 0.003818778 0.003818778 0.003502565 0.004081880
[25] 0.004081880 0.004309115 0.005725876 0.005725876 0.003912183 0.003912183 0.003812053 0.003812053
[33] 0.003446829 0.003446829 0.003985157 0.003985157 0.004181148 0.004181148 0.003772006 0.003642145
[41] 0.003642145 0.003251015 0.003251015 0.003747850 0.003747850 0.003747850 0.003949234 0.003949234
[49] 0.005396190 0.005396190 0.003924456 0.003924456 0.003744212 0.003744212 0.003389661 0.003389661
[57] 0.003744518 0.003744518 0.003744518 0.003962402 0.003962402 0.005306857 0.005306857 0.004970510
[65] 0.004776827 0.004776827 0.004540781 0.004540781 0.004702512 0.004702512 0.004702512 0.004905608
[73] 0.004905608 0.005939234 0.005939234

$df
[1] 69

$residual.scale
[1] 0.015677

diversityDataWB$pred_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$fit
diversityDataWB$se_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$se.fit

plot_pielouJ2_pred <- diversityDataWB %>%

#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_pielouJ2 - 2 * se_pielouJ2, yend = pred_pielouJ2 + 2 * se_pielouJ2, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Evenness (Pielou J2)") #+ scale_y_log10()
plot_pielouJ2_pred

Breakaway richness subplots

plotBreakaway <- diversityDataWB %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  geom_errorbar(aes(ymin = break_lower, ymax = break_upper), width = -.1) + 
  scale_y_log10() +
  ylab("Richness (Breakaway)")
plotBreakaway

#predict(breakLm, se.fit = TRUE)
# doesn't work because built with a different data frame

Why are these not smooth curves?!! What if I redo the model, this time with the same data frame

breakLm2 <- lm(break_estimate ~ log(Size_Class) + I(log(Size_Class) ^2) + lat +  I(lat^2) + depth ,data = diversityDataWB)
breakLm2 %>% summary()

Call:
lm(formula = break_estimate ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityDataWB)

Residuals:
    Min      1Q  Median      3Q     Max 
-2974.5 -1191.2  -151.6   599.9  6768.1 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          7124615.61 3339862.88   2.133   0.0365 *
log(Size_Class)          244.45     202.35   1.208   0.2312  
I(log(Size_Class)^2)     -75.16      37.65  -1.996   0.0498 *
lat                  -370568.38  174035.93  -2.129   0.0368 *
I(lat^2)                4817.28    2265.61   2.126   0.0371 *
depth                    151.10      80.15   1.885   0.0636 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1902 on 69 degrees of freedom
Multiple R-squared:  0.1414,    Adjusted R-squared:  0.0792 
F-statistic: 2.273 on 5 and 69 DF,  p-value: 0.0567

Note the non statistical significance overall

#predict(breakLm2, se.fit = TRUE)
diversityDataWB$pred_break = predict(breakLm2, se.fit = TRUE)$fit
diversityDataWB$se_break = predict(breakLm2, se.fit = TRUE)$se.fit

plot_break_pred <- diversityDataWB %>%

#filter(Depth %in% c("Surface", "Bottom")) %>%
#  filter(Station == 4.3) %>%
  ggplot(aes(x = Size_Class, y = pred_break, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_break - 2 * se_break, yend = pred_break + 2 * se_break, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Richness (Breakaway -- LM)") #+ scale_y_log10()
plot_break_pred

Rebuilding combined products

plotAlphaWB <- plot_grid(plotBreakaway, plotShan, plotPielou2, nrow = 1, labels = LETTERS)
plotAlphaWB

ggsave(here::here("Figures", "BreakawayAlpha.png"), plotAlphaWB, width = 8, height = 4)

Summary table I want both breakaway metrics here

bettaTable <- myBet$table %>% 
  as.data.frame() %>%
  rename(estimate = Estimates,
         `std.error` = `Standard Errors`,
         `p.value`=`p-values`
         ) %>%
  mutate(`statistic` = NA) %>%
  rownames_to_column(var = "term") %>%
  select(term, estimate, std.error, statistic, p.value) %>%
  as_tibble()
bettaTable

alphaSummary2 <- tibble(
  metric = c("Richness (Breakaway -- LM)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
  model = list(breakLm, shanMod, pielouMod2)
)
  
alphaSummary2 <- alphaSummary2 %>%
  mutate(df = map(model, ~broom::tidy(summary(.))))

## Add in willis variables

breakawaySummary <- tibble(
  metric = "Richness (Breakaway -- Betta)",
  model = NULL,
  df = list(bettaTable)
)

alphaSummary2 = bind_rows(breakawaySummary, alphaSummary2)

alphaSummary2 <- alphaSummary2 %>%
  select(-model) %>%
  unnest(df)

alphaSummary2 <- alphaSummary2 %>%
  rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
  mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
         Term = str_replace(Term, "\\^2", "\\^2\\^"),
         Term = str_replace(Term, "depth", "Depth"),
         Term = str_replace(Term, "lat", "Latitude"),
         Term = str_replace(Term, "_", " ")# BOOKMARK!!
         ) %>%
  mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
           reformat_sci()
         ) %>%
  mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
           reformat_sci()
  ) %>%
  mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
  mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
  rename(`Standard\nError` = `Standard Error`) %>%
  identity()



alphaSummary2

alphaTable2 <- alphaSummary2 %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
  align(j = -c(1:2), align = "right")
alphaTable2

Metric	Term	Estimate	Standard Error	T Value	p
Richness (Breakaway – Betta)	Intercept	7.1 x 106	2.4 x 102	NA	< 0.001
	log(Size Class)	1.2 x 102	6.1 x 101	NA	0.058
	log(Size Class)2	-5.0 x 101	1.2 x 101	NA	< 0.001
	Latitude	-3.7 x 105	6.1 x 100	NA	< 0.001
	Latitude2	4.8 x 103	1.6 x 10-1	NA	< 0.001
	Depth	1.5 x 102	1.0 x 101	NA	< 0.001
Richness (Breakaway – LM)	Intercept	7.1 x 106	3.3 x 106	2.13	0.036
	log(Size Class)	2.4 x 102	2.0 x 102	1.21	0.231
	log(Size Class)2	-7.5 x 101	3.8 x 101	-2.00	0.050
	Latitude	-3.7 x 105	1.7 x 105	-2.13	0.037
	Latitude2	4.8 x 103	2.3 x 103	2.13	0.037
	Depth	1.5 x 102	8.0 x 101	1.89	0.064
Diversity (Shannon H)	Intercept	1.4 x 103	7.8 x 102	1.84	0.071
	log(Size Class)	1.9 x 10-1	4.7 x 10-2	4.07	< 0.001
	log(Size Class)2	-3.5 x 10-2	8.7 x 10-3	-4.00	< 0.001
	Latitude	-7.4 x 101	4.0 x 101	-1.83	0.072
	Latitude2	9.6 x 10-1	5.3 x 10-1	1.83	0.072
	Depth	1.5 x 10-2	1.9 x 10-2	0.78	0.436
Evenness (Pielou J)	Intercept	-2.0 x 101	2.8 x 101	-0.72	0.476
	log(Size Class)	-3.3 x 10-3	1.7 x 10-3	-1.97	0.052
	log(Size Class)2	5.8 x 10-4	3.1 x 10-4	1.85	0.068
	Latitude	1.0 x 100	1.4 x 100	0.71	0.477
	Latitude2	-1.3 x 10-2	1.9 x 10-2	-0.71	0.478
	Depth	-2.4 x 10-4	6.6 x 10-4	-0.36	0.718

alphaTable2 %>% save_as_docx(path = here::here("Tables", "alphaTable2.docx"))

myBet$table

And finally predictions from richness, diversity evenness again.

plotAlphaWB_pred <- plot_grid(plot_break_pred,plotShannonH_pred,plot_pielouJ2_pred, nrow = 1, labels = LETTERS)
plotAlphaWB_pred

ggsave(here::here("Figures", "BreakawayAlphaPredictions.png"), plot = plotAlphaWB_pred, width = 8, height = 4)